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Abstract 

We introduce a novel method for numerical spin glass investigations: Simulations of two 
replica at fixed temperature, weighted such that a broad distribution of the Parisi overlap 
parameter q is achieved. Canonical expectation values for the entire g-range (multi-overlap) 
follow by re-weighting. We demonstrate the feasibility of the approach by studying the 3d 
Edwards-Anderson Ising (Jj^ = ±1) spin glass in the broken phase (/3 = 1). For the first 
time it becomes possible to obtain reliable results about spin glass tunneling barriers. In 
addition, as do some earlier numerical studies, our results support that Parisi mean field 
theory is valid down to 3d. 

PACS. 75.40.Mg Numerical simulation studies, 75.50.Lk Spin glasses and other random mag- 
nets 
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One of the questions which ought to be addressed before performing a large scale com- 
puter simulation is "What are suitable weight factors for the problem at hand?" The weight 
factor of canonical Monte Carlo (MC) simulations is exp{—(3E), where E is the energy of 
the configuration to be updated and f3 is the inverse temperature in natural units. The 
Metropolis and other methods generate canonical configurations through a Markov process. 
It has been expert wisdom [|| for quite a while and became widely recognized in recent years 
I, I, 1, 1, 1, 13 that MC simulat ions with a-priori unknown weight factors, like for instance 
the inverse spectral density l/n{E), are also feasible and deserve to be considered. But it is 
not straightforward to design suitable weights. To find a weighting procedure that works in 
practice requires considerable intuitive or other understanding of the underlying physics. 

On the easier side are situations where one has to enhance rare configurations which are 
controlled by some standard thermodynamic observable. They are called static in Ref . . An 
example are temperature-driven first-order phase transitions, controlled by the energy. The 
multicanonical method defines energy dependent weight factors to enhance configurations 
needed to estimate the interfacial tension. The multimagnetical method |Q does the same 
for magnetic field-driven first-order transitions, controlled by the magnetization. Far more 
difficult are problems, called dynamic in Ref. , for which the ergodicity of canonical MC 
simulations breaks down due to energy barriers for which no explicit parameterization in 
terms of a standard thermodynamic variable is known. This is the case for spin glasses and 
other systems with conflicting constraints where the energy barriers are caused by disorder 
and frustration. Below some freezing temperature it becomes extremely difficult to generate 
canonical equilibrium configurations for these systems. Recently, progress has been made by 
exploring p|, ^ |10|, [11| innovative weighting methods for this problem. Along this line our 
paper introduces a novel, efficient approach. 

We focus on the 3d Edwards-Anderson Ising (EAI) spin glass on a simple cubic lattice. 
It is widely considered to be the simplest model to exhibit realistic spin glass behavior and 
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has been the testing ground of Refs. P, p|, |10|, O]. The energy is given by 



E = - J^JikSiSk , (1) 

{ik) 

where the sum is over nearest-neighbour sites. The Ising spins Si and Sk as well as the 
exchange coupling constants Jik take values ±1. A realization is defined by a fixed assign- 
ment of the exchange coupling constants Jik- In our investigation we enforce the constraint 
J2{ik) Jik = by picking half of the Jik at random and assigning them the value +1, whereas 
the others are fixed at —1. 



Early MC simulations of the EAI model, for a concise review see |11|, located the freezing 



temperature at Pc ~ 0.9. Recent, very high statistics canonical simulations |jT2|, |T3| estimate 
Pc = 0.901± 0.034, and considerably improve the evidence in support of a second-order phase 
transition at /3c- The studies [§, ^ ^ focus on first improving the notorious slowing down 
of simulations a.t (3 > (3c- Their main idea is to avoid getting stuck in metastable low-energy 
states by using a Markov process which samples the ordered as well as the disordered regions 
of configuration space in one run. Ref. P] does this by multicanonical |^ re-weighting in 



the energy, whereas Refs. 0, |T0|, |TT[ use and improve the method of enlarged ensembles ^ 
in their simulated tempering version. Refreshing the system in the disordered phase clearly 
benefits the simulations, but the performance has remained below early expectation. It has 
been conjectured that the reason lies in the tree-like structure of the low-energy spin glass 
states, see Ref. for a detailed discussion. 

Studying simulated tempering, Kerler and Rehberg combined two copies (replica) of 
the same realization (defined by its couplings Jik) in one simulation. The purpose was to 
allow for direct evaluation of the Parisi overlap parameter 

^=^1:^hi ■ (2) 
i=i 

Here N denotes the number of spins, the spins s\ = ±1 correspond to the first and the spins 
= ±1 to the second replica. Now, our observation is that one does still control canonical 
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expectation values at temperature (3 ^ when one simulates with a weight function 



w{q) = exp 



PT.Ms]sl + s',sl) + S{q) 



(3) 



This is obvious for S{q) = 0, and a non-trivial S{q) can be mapped onto this situation 
by standard re- weighting ||l|, 0. Of particular interest is to determine S{q) recursively 
such that the histogram H{q) becomes uniform in q and the interpretation of S{q) being 
the microcanonical entropy of the Parisi order parameter. Hence, although an explicit order 
parameter does not exist, an approach very similar to the multimagnetical 0] (which is an 
highly efficient way to sample interface barriers for ferromagnets) exists herewith. In contrast 



to this multi-overlap method, simulated and parallel tempering techniques [T^, |TT| do not 
allow to change barrier heights. 

Our EAI simulations are performed on = lattices at /? = 1, in the interesting 
region well below the freezing temperature. All calculations were done on a cluster of Alpha 
workstations at FSU. We simulated 512 different realizations for L = 4, 6, and 8, and 33 
for L = 12. For all realizations tunneling between the extrema q = ±1 was achieved. Each 
production run of data taking was concluded after at least twenty tunneling event of the 
form 

(g = 0) (g = ±1) and back 

were recorded (for technical reasons the actual numbers were between 20 and 39 per real- 
ization). Table 1 gives an overview of the tunneling performance of our algorithm. Fitting 
the estimates of the mean value r to the form ln(r) = a + z \n{V) gives z = 2.42 ± 0.03. 
Compared with the slowing down of p[ this is an improvement of almost a factor \/V. Still, 
the slowing down is far off from the theoretical optimum z = 1. One reason seems to be 
that we are enforcing the limit g ±1. This limit correlates strongly with ground states. 



which are difficult to reach by local updates, see for instance Being content with a 

smaller region (like the two outmost maxima in the g-distribution) is expected to give fur- 
ther improvements of the tunneling performance. Other data compiled in Table 1 are the 
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encountered minimum, maximum and median tunneling times. We observe that the mean 
values are systematically larger than the median, what means that the tunneling distribution 
has a rather long tail towards large tunneling times. On the other hand, the effect is not 
severely hindering our multi-overlap simulations: For the lattice sizes L = 4 to 8 the worst 
behaved realization took never more than 3% of the entire computer time and for L = 12 
(where we have only 33 realizations) this amount was 12%. 

Initially in each run, a working estimate of the weight function has to be obtained. 



Using a variant of the recursion proposed in 15] this has turned out to be remarkably easy. 



For each case we stopped the recursion of weights after four tunnelings were achieved and 
the used computer time corresponds in good approximation to 4r, with r as given in Table 1. 

The analysis of the thus created data allows us to calculate a number of physically inter- 
esting quantities. In particular accurate determinations of the canonical potential barriers 
in q are, for the first time, possible. Let Pi{q) be the canonical probability densities of q, 
where i = 1, ■■■,n labels the different realizations (additional dependence on lattice size and 
temperature is implicit). We define the corresponding potential barrier by 

-Aq 

Bi= II max [1, Pi{q)/Pi{q + Aq)] , (4) 

where Ag is the stepsize in q. For the double-peak situations of first-order phase transitions 
[|] Eq. (g) becomes Bi = pmax^^min^ where i^™'^'' is the absolute maximum and P™™ is the 
absolute minimum (for ferromagnets at g = 0) of the probability density Pi{q)- Our definition 
generalizes to the situation where several minima and maxima occur due to disorder and 
frustration. When evaluating (^) from numerical data for Pi{q) some care is needed to avoid 
contributions from statistical fluctuations of Pi{q)- 

Graphically, our values for the Bi are presented in Fig. 1. It comes as a surprise that 
the finite-size dependence of the distributions is very weak. Therefore, one may question 
the apparently accepted opinion that these barriers are primarily responsible for the severe 
slowing down of canonical MC simulations with increasing lattice size. To study this issue 
further, we have compiled in Table 2 for each lattice size the following informations about 



our potential barrier results: largest and second largest values -Bmax and B2, median values 
and 70% confidence limits around those, and mean values B with statistical error bars. 
From this table it becomes obvious, why this investigation could not be performed using 
canonical methods to which also enlarged ensembles belong (they enlarge the ensemble but 
still use canonical weights). For these methods the slowing down would be proportional to 
the average barrier height B, which is already large for L = 4, about 18 thousand, and 
increases to about 2.8 million for L = 12. 

Next, the reader may be puzzled by the very large error bars assigned to the mean 
barrier values. Their explanation is: The entire mean value is dominated by the largest 
barrier, which contributes between 70% (L = 4) and, practically, 100% (L = 12), see -Bmax 
in the second column of Table 2. Besides -Bmax? the second largest value B2 is listed in 
the third column. The lesson from these numbers is that very few of the realizations are 
responsible for the collapse of canonical simulation methods. It may be remarked that most 
of these worst case barriers exhibit simple double-peak behavior. An exception is L = 6 
where the distribution yielding -Bmax has two double peaks. To exhibit the difference, the 
inlay of Fig. 1 depicts the right-hand-side of the L = 6 probability densities Pi{q) with 
i = 459 corresponding to the -Bmax and i = 122 to the B2 barrier. For B2 the value of our 
barrier definition (|) agrees with the Pmax/-Pmin value, whereas for -Bmax it is by about a 
factor of two larger. 

Typical configurations, described by the median results of Table 2, have much smaller 
tunneling barriers. They turn out to be quite insensitive to the lattice size, in fact the value 
-Bmed = 12.3 fits into the confidence interval for all simulated lattice sizes. Presumably, there 
is some increase of -Bmed with lattice size, but to trace it we would need to simulate more 
realizations. This result of an almost constant typical tunneling barrier is consistent with 
the fact that our tunneling times are rather far apart from their theoretical optimum: Other 
reasons than overlap barriers have to be responsible. 



Our data are consistent with other numerical evidence 121 in favor of the Parisi mean field 
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scenario being valid down to 3d and against the competing droplet picture, thus indirectly 
supporting Parisi's criticism [|I^ of Ref. |]T6[. The averaged canonical probability densities 



P{q) = [Pi{(l)]scv at the simulation point (3 = 1.0 are shown in Fig. 2. While the peak moves 
with increasing lattice size towards smaller g- values, the value of -P(O) is clearly non-zero and 
shows almost no finite-size dependence. The /3-dependence of P{q) obtained by standard 
reweighting of our time-series data at /9 = 1 is illustrated in Fig. 3 for L = 8. For lower 
temperatures the peak of P{q) becomes more pronounced and moves towards larger g- values. 
The extremal /3-values indicate the inverse temperature range in which reliable results can 
be expected. This range was estimated by measuring the overlap of the reweighted energy 
histogram with the energy histogram at the simulation point f3 = 1.0, individually for each 
of the realizations. With the present statistics the phase transition point should thus be in 
this range, at least up to L = 8. By analyzing the spin glass susceptibility, xsg = N[{q'^)]scv, 
we obtain the best finite-size scaling fit xsg L'^^" at Pc = 0.88 with 7/// = 2.37(4) and 
a goodness-of-fit parameter Q = 0.25. This is corroborated by the curves of the Binder 
parameter, g = (l/2)(3 — [{q*)]a.v/[{q'^)]'iv)y which merge around /? = 0.89. In the low- 
temperature phase (/? > Pc) the curves for different lattice sizes seem to fall on top of each 
other, but our error bars are still too large to draw a firm conclusion from this quantity. 



These results are consistent with the findings of Ref. |[T2[ and could be easily improved by 
redoing the simulations closer to jSc, possibly with more realizations and less statistics per 
realization. For such, or similar, studies the number of L = 12 realizations can be readily 
enhanced by running on a parallel computer like a Cray T3E. Narrowing the g-range will 
allow to simulate lattices of size L = 16 and beyond. 

Finally, we like to mention that our method is particularly well-suited to study the 
infiuence of an interaction term [|1^ 



N 

eY^s]sl = eNq 

i=l 

in the Hamiltonian (|ID: We obtain expectation values for arbitrary e- values. Physically most 
interesting is to combine a non-zero magnetic field with a non-zero e-value. 
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In conclusion, we have demonstrated the feasibihty of using g-dependent (multi-overlap) 
weight factors. Although the tunneling performance is not optimal, the method opens new 
horizons for spin glass simulations. In this paper we succeeded, for the first time, to study q- 
barriers in some details. Using parallel computers and slight modifications of our method (like 
narrowing the g'-range, including a magnetic field, etc.) will allow to extend our investigation 
into various interesting directions, like an improved study of the thermodynamic limit at and 
below the transition point, or e-physics. 
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Tables and Figure Captions 



L 


^min 


^max 


^med 


T 


4 


4.5E02 


6.2E03 


9.9E02 


(1.13±0.03)E03 


6 


4.9E03 


3.1E05 


1.3E04 


(1.88±0.09)E04 


8 


2.4E04 


1.6E06 


1.1E05 


(1.76±0.09)E05 


12 


7.1E05 


1.6E07 


2.7E06 


(4.11±0.65)E06 



Table 1: Overview of the tunneling performance: minimum, maximum, median and mean ± 
error tunneling times. All numbers are in units of sweeps. 



L 


-^max 


B2 


-°med 




-^mcd 


B 


4 


6.56E06 (70%) 


9.11E05 


15.1 


12.4 


9.62 


(1.84± 1.30)E04 


6 


2.76E06 (74%) 


1.44E05 


12.3 


11.1 


10.1 


(7.29 ± 5.42)E03 


8 


1.97E08 (98%) 


1.36E06 


17.7 


15.2 


12.3 


(3.91±3.85)E05 


12 


9.14E07 (100%) 


1.96E03 


35.3 


12.9 


10.7 


(2.77±2.77)E06 



Table 2: Canonical potential barriers: maximum (and its contribution to the mean in %), 
second largest value, upper median confidence limit, median, lower median confidence limit 
(upper and lower limit bound a 70% confidence interval), the mean and its error bar. 



Figure 1: Canonical tunneling barrier distributions at (5 — 1. (The L — 12 barriers are re- 
labelled to fill into the 1-512 range.) The inlay shows the two worst L = 6 realizations. 

Figure 2: Finite-size dependence of the averaged canonical probability densities P{q) at 
(3=1. For L = 8 only every second error bar is shown, and for L = 12 only every 
tenth. 

Figure 3: Temperature dependence of the averaged canonical probability densities P{q) for 
L = 8, obtained by reweighting. Only every second error bar is shown. 
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